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RECORDING MEDIUM, AND SYSTEM ESTIMATION DEVICE 



Technical Field 

5 

The present invention relates to a system estimation method 
and program, a recording medium, and a system estimation device, 
and particularly to a system estimation method and program, a 
recording medium, and a system estimation device, in which the 
10 generation of robustness in state estimation and the optimization 
of a forgetting factor are simultaneously realized by using a fast 
Hop filtering algorithm of a hyper Ha> filter developed on the basis 
of an Hoo evaluation criterion. 

15 Background Art 

In general, system estimation means estimating a parameter 
of a mathematical model (transfer function, impulse response, 
etc. ) of an input/output relation of a system based on input/output 

20 data. Typical application examples include an echo canceller in 
international communication, an automatic equalizer in data 
communication, an echo canceller and sound field reproduction in 
a sound system, active noise control in a vehicle etc. and the 
like. For more information, see non-patent document 1: "DIGITAL 

25 SIGNAL PROCESSING HANDBOOK" 1993, The Institute of Electronics, 
Information and Communication Engineers, and the like. 

(Basic Principle) 

Fig. 8 shows an example of a structural view for system 
30 estimation (unknown system may be expressed by an IIR (Infinite 
Impulse Response) filter) . 

This system includes an unknown system 1 and an adaptive 
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filter 2. The adaptive filter 2 includes an FIR digital filter 
3 and an adaptive algorithm 4. 

Hereinafter, an example of an output error method to identify 
the unknown system 1 will be described. Here, u k denotes an input 
5 of the unknown system 1, d k denotes an output of the system, which 
is a desired signal, and d~ k denotes an output of the filter. 
(Incidentally, means an estimated value and should be placed 
directly above a character, however, it is placed at the upper 
right of the character for input convenience. The same applies 

10 hereinafter . ) . 

Since an impulse response is generally used as a parameter 
of an unknown system, the adaptive filter adjusts a coefficient 
of the FIR digital filter 3 by the adaptive algorithm so as to 
minimize an evaluation error e k = d k - d~ k of the figure. 

15 Besides, conventionally, a Kalman filter based on an update 

expression (Riccati equation) of an error covariance matrix has 
been widely used for the estimation of a parameter (state) of a 
system. The details are disclosed in non-patent document 2: S. 
Haykin: Adaptive filter theory, Prentice-Hall (1996) and the like. 

20 Hereinafter, the basic principle of the Kalman filter will 

be described. 

A minimum variance estimate x~ k | k of a state x k of a linear 
system expressed in a state space model as indicated by the 
following expression: 
25 x k+1 = p" 1/2 x k , y k = H k x k + v k (1) 

is obtained by using an error covariance matrix 2)~k|k-i of the state 
as follows. 
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*k\k = ^k\k-i + K k (y k -H k x k \ k .i) 



i 




(2) 
(3) 



(4) 



where , 



x 0 |_i = 0, i7 0 |_i = e 0 I, e 0 > 0 



(5) 



x k : State vector or simply a state; unknown and this is an object 
of estimation. 

y k : Observation signal; input of a filter and known. 
5 H k : Observation matrix; known. 
v k : Observation noise; unknown. 

p: Forgetting factor; generally determined by trial and error. 
K k : Filter gain; obtained from matrix 2)~ k | k -i. 

2~ k | k : Corresponds to the covariance matrix of an error of x~ k | k ; 
10 obtained by a Riccati equation. 

2)%+i|k: Corresponds to the covariance matrix of an error of x" k+ i| k ; 
obtained by the Riccati equation. 

2~i| 0 : Corresponds to the covariance matrix in an initial state; 
although originally unknown, e 0 I is used for convenience. 

15 The present inventor has already proposed a system 

identification algorithm by a fast Ha filter (see patent document 
1). This is such that an Hoo evaluation criterion is newly 
determined for system identification, and a fast algorithm for 
the hyper H« filter based thereon is developed, while a fast 

20 time-varying system identification method based on this fast H* 
filtering algorithm is proposed. The fast Hoo filtering algorithm 
can track a time-varying system which changes rapidly with a 
computational complexity of O (N) per unit -time step. It matches 
perfectly with a fast Kalman filtering algorithm at the limit of 
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the upper limit value. By the system identification as stated 
above, it is possible to realize the fast real-time identification 
and estimation of the time- invariant and time-varying systems. 
Incidentally, with respect to methods normally known in the 
5 field of the system estimation, see, for example, non-patent 
documents 2 and 3. 

(Applied Example to Echo Canceller) 

In a long distance telephone circuit such as an international 

10 telephone, a four-wire circuit is used from the reason of signal 
amplification and the like. On the other hand, since a 
subscriber's circuit has a relatively short distance, a two-wire 
circuit is used. 

Fig. 9 is an explanatory view concerning a communication 

15 system and an echo. A hybrid transformer as shown in the figure 
is introduced at a connection part between the two -wire circuit 
and the four-wire circuit, and impedance matching is performed. 
When the impedance matching is complete, a signal (sound) from 
a speaker B reaches only a speaker A. However, in general, it is 

20 difficult to realize the complete matching, and there occurs a 
phenomenon in which part of the received signal leaks to the 
four-wire circuit, and returns to the receiver (speaker A) after 
being amplified. This is an echo (echo). As a transmission 
distance becomes long (as a delay time becomes long) , the influence 

25 of the echo becomes large, and the quality of a telephone call 
is remarkably deteriorated (in the pulse transmission, even in 
the case of short distance, the echo has a large influence on the 
deterioration of a telephone call). 

Fig. 10 is a principle view of an echo canceller. 

30 Then, as shown in the figure, the echo canceller (echo 

canceller) is introduced, an impulse response of an echo path is 
successively estimated by using a received signal which can be 
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directly observed and an echo, and a pseudo-echo obtained by using 
it is subtracted from the actual echo to cancel the echo and to 
remove it. 

The estimation of the impulse response of the echo path is 
5 performed so that the mean square error of a residual echo e k becomes 
minimum. At this time, elements to interfere with the estimation 
of the echo path are circuit noise and a signal (sound) from the 
speaker A. In general, when two speakers simultaneously start to 
speak (double talk), the estimation of the impulse response is 
10 suspended. Besides, since the impulse response length of the 
hybrid transformer is about 50 [ms], when the sampling period is 
made 125 [\is] , the order of the impulse response of the echo path 
becomes actually about 400. 
Non-patent document 1 
15 "DIGITAL SIGNAL PROCESSING HANDBOOK" 1993 The Institute of 

Electronics, Information and Communication Engineers 
Non -patent document 2 

S. Haykin: Adaptive filter theory. Prentice-Hall (1996) 
Non-patent document 3 
20 B. Hassibi, A. H. Sayed, and T. Kailath: 

"Indefinite-Quadratic Estimation and Control", SIAM (1996) 
Patent document 1 
JP-A-2002-135171 

25 Disclosure of the Invention 

However, in the conventional Kalman filter including the 
forgetting factor p as in the expressions (1) to (5), the value 
of the forgetting factor p must be determined by trial and error 
30 and a very long time has been required. Further, there has been 
no means for judging whether the determined value of the forgetting 
factor p is an optimal value. 
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Besides, with respect to the error covariance matrix used 
in the Kalman filter, it is known that a quadratic form to an 
arbitrary vector, which is originally not zero, is always positive 
(hereinafter referred to as "positive definite") , however, in the 
5 case where calculation is performed by a computer at single 
precision, the quadratic form becomes negative (hereinafter 
referred to as "negative definite"), and becomes numerically 
unstable. Besides, since the amount of calculation is 0(N 2 ) (or 
0(N 3 )), in the case where the dimension N of the state vector x k 

10 is large, the number of times of arithmetic operation per time 
step is rapidly increased, and it has not been suitable for a 
real-time processing. 

In view of the above, the present invention has an object 
to establish an estimation method which can theoretically 

15 optimally determine a forgetting factor, and to develop an 
estimation algorithm and a fast algorithm which are numerically 
stable. Besides, the invention has an object to provide a system 
estimation method which can be applied to an echo canceller in 
a communication system or a sound system, sound field reproduction, 

20 noise control and the like. 

In order to solve the problem, according to the invention, 

a newly devised Hoo optimization method is used to derive a state 
estimation algorithm in which a forgetting factor can be optimally 
determined. Further, instead of an error covariance matrix which 

25 should always have the positive definite, its factor matrix is 
updated, so that an estimation algorithm and a fast algorithm, 
which are numerically stable, are developed. 

According to first solving means of the invention, 

a system estimation method and program and a computer 

30 readable recording medium recording the program are for making 
state estimation robust and optimizing a forgetting factor p 
simultaneously in an estimation algorithm, in which 
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for a state space model expressed by following expressions; 
x k+1 = F k x k + G k w k 
y k = H k x k + v k 
z k = H k x k 
5 here , 

x k : a state vector or simply a state, 

w k : a system noise, 

v k : an observation noise, 

y k : an observation signal, 
10 z k : an output signal, 

F k : dynamics of a system, and 

G k : a drive matrix, 

a maximum energy gain to a filter error from a disturbance 

weighted with the forgetting factor p as an evaluation criterion 
15 is suppressed to be smaller than a term corresponding to a 

previously given upper limit value y f , and 

the system estimation method, the system estimation program 

for causing a computer to execute respective steps, and the 

computer readable recording medium recording the program, 
20 includes 

a step at which a processing section inputs the upper limit 

value Yfr the observation signal y k as an input of a filter and 

a value including an observation matrix H k from a storage section 

or an input section, 
25 a step at which the processing section determines the 

forgetting factor p relevant to the state space model in accordance 

with the upper limit value y f , 

a step at which the processing section reads out an initial 

value or a value including the observation matrix H k at a time from 
30 the storage section and uses the forgetting factor p to execute 

a hyper H« filter expressed by a following expression: 

x~ k | k = F k -ixVi| k _i + K s , k (y k - HjcFjc.iX^^!)^!) 
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here , 

x~ k | k ; an estimated value of a state x k at a time k using observation 
signals y 0 to y k/ 
F k : dynamics of the system, and 
5 K s , k ; a filter gain, 

a step at which the processing section stores an obtained 
value relating to the hyper H„ filter into the storage section, 

a step at which the processing section calculates an 
existence condition based on the upper limit value Yf and th © 
10 forgetting factor p by the observation matrix Hi and the filter 
gain K s#i , and 

a step at which the processing section sets the upper limit 

value to be small within a range where the existence condition 

is satisfied at each time and stores the value into the storage 
15 section by decreasing the upper limit value Yf and repeating the 

step of executing the hyper filter. 

Besides, according to second solving means of the invention, 
a system estimation device is for making state estimation 

robust and optimizing a forgetting factor p simultaneously in an 
20 estimation algorithm, in which 

for a state space model expressed by following expressions: 

x k+1 = F k x k + G k w k 

y k = H k x k + v k 

Zk = H k x k 
25 here, 

x k ; a state vector or simply a state, 

w k : a system noise, 

v k : an observation noise, 

y k : an observation signal, 
30 z k : an output signal, 

F k : dynamics of a system, and 

G k : a drive matrix. 
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a maximum energy gain to a filter error from a disturbance 
weighted with the forgetting factor p as an evaluation criterion 
is suppressed to be smaller than a term corresponding to a 
previously given upper limit value Yf* 
5 the system estimation device includes: 

a processing section to execute the estimation algorithm; 

and 

a storage section to which reading and/or writing is 
performed by the processing section and which stores respective 
10 observed values, set values, and estimated values relevant to the 
state space model, 

the processing section inputs the upper limit value Yf , the 
observation signal y k as an input of a filter and a value including 
an observation matrix H k from the storage section or an input 
15 section, 

the processing section determines the forgetting factor p 
relevant to the state space model in accordance with the upper 
limit value y f , 

the processing section reads out an initial value or a value 
20 including the observation matrix H k at a time from the storage 
section and uses the forgetting factor p to execute a hyper H«> filter 
expressed by a following expression: 

x~ k | k = F k -ixVi| k -i + K s , k (y k - HjcFk-ixViM) 
here, 

25 x~ k | k ; an estimated value of a state x k at a time k using observation 
signals y 0 to y k , 
F k : dynamics of the system, and 
K Sfk ; a filter gain, 

the processing section stores an obtained value relating 
30 to the hyper H«> filter into the storage section, 

the processing section calculates an existence condition 
based on the upper limit value Yf and the forgetting factor p by 
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the observation matrix Hi and the filter gain K s#1 , and 

the processing part sets the upper limit value to be small 
within a range where the existence condition is satisfied at each 
time and stores the value into the storage section by decreasing 
5 the upper limit value Yf and repeating the step of executing the 
hyper filter. 

According to the estimation method of the invention, the 
forgetting factor can be optimally determined, and the algorithm 
can stably operate even in the case of single precision, and 
10 accordingly, high performance can be realized at low cost. In 
general, in a normal civil communication equipment, calculation 
is often performed at single precision in view of cost and speed. 
Thus, as the practical state estimation algorithm, the invention 
would have effects in various industrial fields. 

15 

Brief Description of the Drawings 

Fig. 1 is a structural view of hardware of an embodiment. 
Fig. 2 is a flowchart concerning the generation of robustness 
20 of an H«> filter and the optimization of a forgetting factor p. 

Fig. 3 is a flowchart of an algorithm of the H«> filter (S105) 
in Fig. 2. 

Fig. 4 is an explanatory view of a square root array algorithm 
of Theorem 2. 

25 Fig. 5 is a flowchart of a fast algorithm of Theorem 3, which 

is numerically stable. 

Fig. 6 is a view showing values of an impulse response 
{hjji-o 23 * 

Fig. 7 shows an estimation result of the impulse response 
30 by the fast algorithm of Theorem 3, which is numerically stable. 
Fig. 8 is a structural view for system estimation. 
Fig. 9 is an explanatory view of a communication system and 
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an echo. 

Fig. 10 is a principle view of an echo canceller. 
DESCRIPTION OF THE PREFERED EMBODIMENTS 

5 

Hereinafter, embodiments of the invention will be described. 
1. Explanation of Symbols 

10 First, main symbols used in the embodiments of the invention 

and whether they are known or unknown will be described. 

x k : State vector or simply a state; unknown and this is an object 

of the estimation. 

x 0 : Initial state; unknown. 
15 w k : System noise; unknown. 

v k : Observation noise; unknown. 

y k : Observation signal; input of a filter and known. 
z k : Output signal; unknown. 
F k : Dynamics of a system; known. 
20 G k : Drive matrix; known at the time of execution. 
H k : Observation matrix; known. 

x~ k | k : Estimated value of a state x k at a time k, using observation 
signals y 0 to y k ; given by a filter equation. 

x" k+ i| k : Estimated value of a state x k+ i at a time k+1 using the 
25 observation symbols y 0 to y k ; given by the filter equation. 

x~ 0 |o: Initial estimated value of a state; originally unknown, 
however, 0 is used for convenience. 

2~k|k: Corresponds to a covariance matrix of an error of x~ k | k ; given 
by a Riccati equation. 

30 2~ k+ i| k : Corresponds to a covariance matrix of an error of x~ k +i| k ; 
given by the Riccati equation. 

2~i| 0 : Corresponds to a covariance matrix in an initial state; 
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originally unknown, however, e 0 I is used for convenience. 
K s , k : Filter gain; obtained from a matrix 2~k| k -i. 
p; Forgetting factor; in the case of Theorems 1 to 3, when Yf is 
determined, it is automatically determined by p = 1 - %(Yf)- 
5 e f<i : Filter error 

R e ,k: Auxiliary variable 

Incidentally, n ~" and "v" placed above the symbol mean 
estimated values . Besides, n - n , "U w and the like are symbols 

added for convenience. Although these symbols are placed at the 

10 upper right of characters for input convenience, as indicated in 
mathematical expressions, they are the same as those placed 
directly above the characters. Besides, x, w and the like are 
vectors , H, G, K, R, J and the like are matrixes and should be 
expressed by thick letters as indicated in the mathematical 

15 expressions, however, they are expressed in normal letters for 
input convenience. 

2. Hardware and Program for System Estimation 

20 The system estimation or the system estimation device and 

system can be provided by a system estimation program for causing 
a computer to execute respective procedures, a computer readable 
recording medium recording the system estimation program, a 
program product including the system estimation program and 

25 capable of being loaded into an internal memory of a computer, 
a computer, such as a server, including the program, and the like. 
Fig. 1 is a structural view of hardware of this embodiment. 
This hardware includes a processing section 101 which is 
a central processing unit (CPU), an input section 102, an output 

30 section 103, a display section 104 and a storage section 105. 
Besides, the processing section 101, the input section 102, the 
output section 103, the display section 104 and the storage section 
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105 are connected by suitable connection means such as a star or 
a bus. Known data indicated in "1. Explanation of Symbols" and 
subjected to the system estimation are stored in the storage 
section 105 as the need arises. Besides, unknown and known data, 
5 calculated data relating to the hyper H» filter, and other data 
are written and/or read by the processing section 101 as the need 
arises. 

3. Hyper H w Filter by Which Forgetting Factor Can Be Optimally 
10 Determined 

(Theorem 1) 

Consideration is given to a state space model as indicated 
by following expressions. 

= F k x k + G k w k) w k , x k 6 Tl N (6) 
y k = H k x k + v k , y*,v*€ft (7) 
z k = H k x k , z k eTZ, H k eTZ ]xN , A: = 0, 1, . . . , L (8) 



15 



An Hoo evaluation criterion as indicated by the following 
expression is proposed for the state space model as described 
above . 

i \m\Vp 

1=0 . 9 



SU P, " i ... k < v 

(9) 



x *' {WiUvi) ll*o - *oi-ill> . + £ INI 2 + £ NIVp 

t=0 t=0 



20 A state estimated value x~ k | k (or an output estimated value 

z v k|k) to satisfy this Hoo evaluation criterion is given by a hyper 
Hoo filter of level Yf: 
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h\k = H k x k \ k 

x k \ k = J 7 * *-i«*:_i|jb_i + K $tk (y k - H k F k ^x k _ l]k _i) 
K s , k = E^HliHkS^Hl + p)- 1 
%k\k = %k\k-i - ^k\k-\C\R~lCkE k \k-\ 
%k+l\k = ( F k£k\kFl) lp 



where , 



e /,« = ~ H&u #0|0 = £l|0 = £o 



R e ,k = Rk + C k 2k\k-\Cl, Rk = 
0<p = l-x(7/)<l, 7/>l 



0 -r- 2 



0 

■PI) 



C k = 



H k 



(10) 

(11) 
(12) 



} (13) 



(14) 
(15) 



Incidentally, expression (11) denotes a filter equation, 
expression (12) denotes a filter gain, and expression (13) denotes 
5 a Riccati equation. 

Besides, a drive matrix G k is generated as follows, 

GtG T = xMp^Fj (16) 

Besides, in order to improve the tracking capacity of the 
foregoing H„ filter, the upper limit value Yf is set to be as small 
10 as possible so as to satisfy the following existence condition. 

E~l = z-^ + I^iIhJH^O, t = 0,...,fc (17) 

ft 

Where, x(Yf) is a monotonically damping function of y t , which 
satisfies x(!) = 1 and x(°°) = 0- 

The feature of Theorem 1 is that the generation of robustness 
15 in the state estimation and the optimization of the forgetting 
factor p are simultaneously performed. 

Fig. 2 shows a flowchart concerning the generation of 
robustness of the H«> filter and the optimization of the forgetting 



14 



P0737PC_Art 34_Amendment050609 



factor p. Here, 

block "EXC > 0": an existence condition of the H«> filter, and 
Ay: a positive real number. 

First, the processing section 101 reads out or inputs the 
5 upper limit value y t from the storage section 105 or the input 
section 102 (S101). In this example, Yf >> 1 is given. The 
processing section 101 determines the forgetting factor p by- 
expression (15) (S103). Thereafter, the processing section 101 
executes the hyper filter of expression (10) to expression (13) 

10 based on the forgetting factor p (S105). The processing section 
101 calculates expression (17) (or the right side (this is made 
EXC) of after-mentioned expression (18)) (S107), and when the 
existence condition is satisfied at all times (S109), Yf is 
decreased by Ay, and the same processing is repeated (Sill). On 

15 the other hand, when the existence condition is not satisfied at 
a certain Yf (S109) , what is obtained by adding Ay to the Yf is made 
the optimal value Yf op of Yf , and is outputted to the output section 
103 and/or stored into the storage section 105 (S113). 
Incidentally, in this example, although Ay is added, a previously 

20 set value other than that may be added. This optimization process 
is called ay-iteration. Incidentally, the processing section 101 
may store a suitable intermediate value and a final value obtained 
at respective steps, such as the H» filter calculation step S105 
and the existence condition calculation step S107 , into the storage 

25 section 105 as the need arises, and may read them from the storage 
section 105. 

When the hyper Hoc filter satisfies the existence condition, 
the inequality of expression (9) is always satisfied. Thus, in 
the case where the disturbance energy of the denominator of 
30 expression (9) is limited, the total sum of the square estimated 
error of the numerator of expression (9) becomes bounded, and the 
estimated error after a certain time becomes 0. This means that 
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when Yf can be made smaller, the estimated value x~ k | k can quickly 
follow the change of the state x k . 

Here, attention should be given to the fact that the 
algorithm of the hyper H« filter of Theorem 1 is different from 
5 that of the normal H«, filter. Besides, when Yf -* 00 , then p = 1 
and G k ■ 0 , and the algorithm of the H* filter of Theorem 1 coincides 
with the algorithm of the Kalman filter. 

Fig. 3 is a flowchart of the algorithm of the (hyper) filter 
(S105) in Fig. 2. 

10 The hyper H<x> filtering algorithm can be summarized as 

follows . 

[Step S201] The processing section 101 reads out the initial 
condition of a recursive expression from the storage section 105 
15 or inputs the initial condition from the input section 102, and 
determines it as indicated in the figure. Incidentally, L denotes 
a previously fixed maximum data number. 

[Step S203] The processing section 101 compares the time k with 
the maximum data number L . When the time k is larger than the 

20 maximum data number, the processing section 101 ends the processing, 
and when not larger, advance is made to a next step. (If 
unnecessary, the conditional sentence can be removed. 
Alternatively, restart may be made as the need arises.) 
[Step S205] The processing section 101 calculates a filter gain 

25 K s#k by using expression (12). 

[Step S207] The processing section 101 updates the filter 
equation of the hyper H<x> filter of expression (11). 
[Step S209] The processing section 101 calculates terms 2~ k | k , 
X*k+i|k corresponding to the covariance matrix of an error by using 

30 the Riccati equation of expression (13). 

[Step S211] The time k is made to advance (k = k + 1) , return is 
made to step S203, and continuation is made as long as data exists. 
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Incidentally, the processing section 101 may store a 
suitable intermediate value, a final value, a value of the 
existence condition and the like obtained at the respective steps, 
such as the H*, filter calculation steps S205 to S209, into the 
storage section 105 as the need arises, or may read them from the 
storage section 105. 

(Scalar Existence Condition) 

The amount of calculation 0(N 2 ) was necessary for the 
judgment of the existence condition of expression (17) . However, 
when the following condition is used, the existence of the H* filter 
of Theorem 1, that is, expression (9) can be verified by the amount 
of calculation O(N) . 

Corollary 1: Scalar existence condition 

When the following existence condition is used, the 
existence of the hyper H<x> filter can be judged by the amount of 
calculation O(N) . 

+ P7/>0, 1 = 0, — ,* (1 8 ) 

Here, 

6 = *= 1 P -h!k,S p = l ~^f) (19) 

Where, K s<1 denotes the filter gain obtained in expression (12). 

(Proof) 

Hereinafter, the proof of the system 1 will be described. 
When a characteristic equation 
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eM -H k E k{k _ x H T k A - {-py} + HtZ^BZ) 

= A 2 - {2H k E k ^ x H T k + pg)X - p^j + pQHt&w^Bt! = 0 

of a 2 x 2 matrix R e ,k is solved, an eigenvalue X± of R e ,k is obtained 
as follows. 

, # ± - Ap e H k E k \ k . x H T k + 4/^7? 
A, 



Where , # = 2H k E k ^ x Hl + PQ,Q=l-i} 



If 



5 one of two eigenvalues of the matrix R e ,k becomes positive, the 
other becomes negative, and the matrixes R k and R e ,k have the same 
inertia. By this, when 

H k E k \ k _ i H k r = H k k k = 

is used, the existence condition of expression (18) is obtained. 
10 Here, the amount of calculation of H k K s , k is O(N) . 

4. State Estimation Algorithm Which Is Numerically Stable 



Since the foregoing hyper H„ filter updates 2~k|k-i G R nxn , 
15 the amount of calculation per unit time step becomes 0(N 2 ), that 
is, an arithmetic operation proportional to N 2 becomes necessary. 
Here, N denotes the dimension of the state vector x k . Thus, as 
the dimension of x k is increased, the calculation time required 
for execution of this filter is rapidly increased. Besides, 
20 although the error covariance matrix 2"k|k-i must always have the 
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positive definite from its property, there is a case where it has 
numerically the negative definite. Especially , in the case where 
calculation is made at single precision, this tendency becomes 
remarkable. At this time # it is known that the filter becomes 
unstable. Thus, in order to put the algorithm to practical use 
and to reduce the cost, the development of the state estimation 
algorithm which can be operated even at single precision (example: 
32 bit) is desired. 

Then, next, attention is paid to 

10 R k =R 1/2 k JiR^ /2 k, 

p _pl/2 t pTi/2 
Ke,k-K e,k*JlK- e ,k* 

o k|k-l=o k|k-lS - k |k-l 

and an Hoo filter (square root array algorithm) of Theorem 1, which 
15 is numerically stabilized, is indicated in Theorem 2. Here, 
although it is assumed that F k = I is established for simplification, 
it can be obtained in the same way also in the case of F k * I. 
Hereinafter, the hyper Hoo filter to realize the state estimation 
algorithm which is numerically stable will be indicted. 



20 



(Theorem 2) 
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*k|fc = *t-i|t-i + K fJe (y k - H k x k .i\ k -i) 
K.t-Kki:, 1W1, 1) , K^p^K&lJ?)!^ 



Rl 



„ i 



0 ^t|t-l J 



e(*) = 



1 

^4 — 


0 







(20) 
(21) 

(22) 



Where, 

Rt^RlJiRl Rl = 
Re* = Rk + Ct^in-iC^, Ct = 



r /9i 0 
0 P*7/. 



1 0 
0 -1 



, £fc|*-i - El\k-iEk\k-i 



l A 



(23) 



0(k) denotes a J-unitary matrix, that is, satisfies 
0(k)J0H(k) T = J, J = (J x 0 I), I denotes a unit matrix, 
K k (:,l) denotes a column vector of a first column of the 
matrix K k . 

Incidentally, in expressions (21) and (22), Ji" 1 and J x 
can be deleted. 

Fig. 4 is an explanatory view of the square root array 
5 algorithm of Theorem 2. This calculation algorithm can be used 
in the calculation (S105) of the H*, filter in the flowchart of 
Theorem 1 shown in Fig. 2. 

In this estimation algorithm, instead of obtaining 2~k|k-i 
by a Riccati type update expression, its factor matrix 2" 1/2 k|k-i 
10 G R NXN (square root matrix of 2~k|k-i) is obtained by the update 
expression based on the J-unitary transformation. From a 1-1 
block matrix and a 2-1 block matrix generated at this time, the 
filter gain K s , k is obtained as shown in the figure. Thus, 2~k|k-i 
= 2~ 1/2 k!k-i S^^kik-i > 0 is established, the positive definite 
15 property of 2~k|k-i is ensured, and it can be numerically stabilized. 
Incidentally, a computational complexity of the H«> filter of 
Theorem 2 per unit step remains 0(N 2 ) . 

Incidentally, in Fig. 4, Ji" 1 can be deleted. 
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10 



First, the processing section 101 reads out terms contained 
in the respective elements of the left -side equations of expression 
(22) from the storage section 105 or obtains them from the internal 
memory or the like, and executes the J-unitary transformation 
(S301) . The processing section 101 calculates system gains K k and 
K s , k from the elements of the right-side equations of the obtained 
expression (22) based on expression (21) (S303, S305). The 
processing section 101 calculates the state estimated value x~ k | k 
based on expression (20) (S307). 

5. Numerically Stable Fast Algorithm for State Estimation 



As described above, a computational complexity of the Hoo 
filter of Theorem 2 per unit step remains 0(N 2 ). Then, as a 
15 countermeasure for the complexity, by using that when H k = H k+ iV, 
H k ■ [u(k), • • •, u(0), 0, • • 0], a covariance matrix 2 k+ i k | k of 
one step prediction error of x k = [x T k , 0 T ] T satisfies 



0 



(24) 



consideration is given to updating L k (that is, L~ k ) with a low 
20 dimension instead of 2 k+1 | k . Here, when attention is paid to R r , k 
| ■ R 1/2 r ,k SR~ i/2 r , k , next Theorem 3 can be obtained. 

(Theorem 3) 
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e(jt) (63) 



**|* = + K t jk{y k - H*x*_i|fc-i) (61) 

K,t = JTfcO, 1)//^(1, 1) , iC fc = pl^JCjjHt (62) 

here, 0(k) denotes an arbitrary J -unitary matrix, and C k = 
Ck4iW is established, where 



0 p 



4]' Ji= [5 -°i 



(23) 



Incidentally, the proof of Theorem 3 will be described later. 
The above expression can be arranged with respect to K k 
instead of K" k (= P" 1/2 K k ) . 

Further, when the following J -unitary matrix 

®(fc)=cy?l © -4 )2 (*)(<t^ ©-<b 

is used, a fast state estimation algorithm of Theorem 4 can be 
obtained. Where, *P denotes a shift matrix. 



(Theorem 4) 
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0 


0 







o 



#r,*+l = -Rr.it - ^JfeCfc+l^e.Jk^+l^'ik 



(25) 
(26) 

(27) 

(28) 

(29) 
(30) 



Where, 



Hk+\ 



fie,l = Hi + C 1 i?i|oC 1 , flj = 



1 0 
0 0 
0 1 



H k+ x = [u k+l u(k + l- N)} = [u(k + 1) u*l, Hi = [u(l), 0, . . . , 0] 

, i?i|o = diag{p 2 )/0 3 ,...,p^ }) p= i_ x(7/) 

, Ko = 0, i 0)0 = x 0 , K t = p" (31) 



P 0 
0 -pyj 



-1 0 



diag[ ] denotes a diagonal matrix, and R e .k+i(l,l) denotes a 1-1 
element of the matrix R e ,k+i. Besides, the above expression can 
be arranged with respect to instead of K" k . 

In the fast algorithm, since the filter gain K s , k is obtained 
by the update of L~ k e r (n+1 ' x2 in the following factoring 



(32) 



0(N+1) is sufficient for the amount of calculation per unit step. 
Here, attention should be paid to the following expression. 



Kk+i 




0 


0 







Fig. 5 is one example of a flowchart of a numerically stable 
fast algorithm of Theorem 34. The fast algorithm is incorporated 
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in the calculation step (S105) of the H«> filter of Fig. 2, and 
is optimized by the y-iteration. Thus, during a period in which 
the existence condition is satisfied, Yf is gradually decreased, 
however, at the time point when it comes to be unsatisfied, Yi is 
increased as indicated in the figure. 

The Hoo filtering algorithm can be summarized as follows. 

[Step S401] The processing section 101 determines an initial 
condition of the recursive expression as indicated in the figure. 
Incidentally, L denotes a maximum data number. 
[Step S403] The processing section 101 compares the time k with 
the maximum data number L. When the time k is larger than the 
maximum data number, the processing section 101 ends the processing, 
and when not larger, advance is made to a next step. (When 
unnecessary, the conditional sentence can be removed. 
Alternatively, restart is made.) 

[Step S405] The processing section 101 recursively calculates a 
term K k+1 corresponding to a filter gain by using expressions (27) 
and (31) . 

[Step S406] The processing section 101 recursively calculates 
Re,k+i by using expression (29). 

[Step S407] The processing section 101 further calculates K s ,k by 
using expressions (26) and (31). 

[Step S409] The processing section 101 judges the existence 
condition EXC > 0 here, and when the existence condition is 
satisfied, advance is made to step S411. 

[Step S413] On the other hand, when the existence condition is 
not satisfied at step S409, the processing section 101 increases 
Yf, and return is made to step S401. 

[Step S411] The processing section 101 updates the filter 
equation of the Hoo filter of expression (25). 

[Step S415] The processing section 101 recursively calculates 
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Rr,k+i by using expression (30). Besides, the processing section 
101 recursively calculates L\ +i by using expressions ( 28 ) and ( 31 ) . 
[Step S419] The processing section 101 advances the time k (k = 
k+1), returns to step S403, and continues as long as data exists. 

5 

Incidentally, the processing section 101 may store a 
suitable intermediate value and a final value obtained at the 
respective steps, such as the H* filter calculation steps S405 
to S415 and the calculation step S409 of the existence condition, 
10 into the storage section 105 as the need arises, and may read them 
from the storage section 105. 

6. Echo Canceller 

15 Next, a mathematical model of an echo canceling problem is 

generated. 

First, when consideration is given to the fact that a 
received signal {u k } becomes an input signal to an echo path, by 
a (time- varying) impulse response {h ± [k]> of the echo path, an 
20 observed value {y k } of an echo {d k } is expressed by the following 
expression. 

N-l 

Vk = d k + v k = ^,hi[k]u k -i + v ki A; = 0,1,2,"- (33) 

t=0 

Here, u k and y k respectively denote the received signal and the 
echo at a time t k (= kT; T is a sampling period) , v k denotes circuit 
25 noise having a mean value of 0 at the time t k , h ± [k] , i = 0, • • 
N-l denotes a time-varying impulse response, and the tap number 
N thereof is known. At this time, when estimated values {h~i[k]> 
of the impulse response are obtained every moment, a pseudo-echo 
as indicated below can be generated by using that. 
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dk = £M*R-i, = 0,1,2,.. (34) 

When this is subtracted from the echo (y k - d~ k « 0) , the echo can 
be cancelled. Where, it is assumed that if k-i < 0, then u k -i = 
0. 

5 From the above, the problem can be reduced to the problem 

of successively estimating the impulse response {hi[k] } of the echo 
path from the received signal {u k } and the echo {y k } which can be 
directly observed . 

In general, in order to apply the H«, filter to the echo 
10 canceller, first, expression (32) must be expressed by a state 
space model including a state equation and an observation equation. 
Then, since the problem is to estimate the impulse response {h ± [k] }, 
when {h ± [k]} is made a state variable x k , and a variance of about 
w k is allowed, the following state space model can be established 
15 for the echo path. 

Xfc+i = x k + G k w k , x k ,w k eH N (35) 

y k = H k x k + v k > VkiVkZn (36) 

Zk = H h x k , z k e 71, H k e n lxN (37) 

Where , 

x k = [h<>[k]r-->h N - X [k}] T , w k = [w k {l)r--iMN)) T 

H k = [ujt,---,Ujfe_tf + i] 

The hyper and fast H«> filtering algorithms to the state space 
model as stated above is as described before. Besides, at the 
estimation of the impulse response, when the generation of a 
20 transmission signal is detected, the estimation is generally 
suspended during that. 

7. Evaluation to Impulse Response 
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(Confirmation of Operation) 

With respect to the case where the impulse response of the 
echo pulse is temporally invariable (hi[k] = hi) , and the tap number 
5 N is 48, the operation of the fast algorithm is confirmed by using 
a simulation. 

47 

Vk = Y, h i u k-i + V k 

«=° (38) 

Incidentally, Fig. 6 is a view showing values of the impulse 

response {hi} here. 
10 Here, the value shown in the figure are used for the impulse 

response {hi} iB0 23 , and the other {hi}i„ 2 4 47 is made 0. Besides, it 

is assumed that u k is stationary Gaussian white noise having a mean 

value of 0 and variance a v 2 = 1.0 x 10" 6 , and the sampling period 

T is made 1.0 for convenience. 
15 Besides, the received signal {u k } is approximated by a 

secondary AR model as follows. 

u k = aiu k -i + <x 2 u k - 2 + w k ' (39) 

Where, a x = 0.7 and ct 2 » 0.1 are assumed, and w k ' denotes stationary- 
Gaussian white noise having a means value of 0 and variance o w . 2 
20 = 0.04. 

(Estimation Result of Impulse Response) 

Fig. 7 shows an estimation result of the impulse response 
| by the numerically stable fast algorithm of Theorem 43. Here, the 
25 vertical axis of Fig. 7(b) indicates 
^{W 7 (hi - x~ k (i + l)) 2 }. 

By this, it is understood that the estimation can be 
excellently performed by the fast algorithm. Where, 

p = 1 - x(Yf)* X(Yf) - Yf~ 2 ' x~ 0 |o = 0 and 2~o|o = 201 were assumed, 
30 and the calculation was performed at double precision. Besides, 
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while the existence condition is confirmed, Yf = 5 - 5 was set 
8. Proof of Theorem 

8-1. Proof of Theorem 2 

When the following expression: 



R e* 



T 

RJ 0 



J k\k- 

0 



(40) 



10 



is established, following expressions are obtained by comparing 
the respective terms of 2 x 2 block matrixes of both sides. 

Re,k = Rk + C k S k \k-iC\ (41) 
= i^iA-iCjf (42) 
+ p-'KkR^Kl = p- 1 !?*,^ (43) 

This is coincident with the Riccati equation of expression 
(13) at Fit = I of Theorem 1. Where, 



J = (Jiei), j 1 = 



1 0 
0 -1 



c k = 



H k 



(44) 



On the other hand, when AJA T = BJB T is established, B can 
15 be expressed asB = A0(k) by using the J-unitary matrix 0(K) . Thus, 
from expression (40), the Riccati equation of Theorem 1 is 
equivalent to the following expression. 



R} 



o 



e(k) 



(45) 



Incidentally, in expressions (40) and (45), Ji" 1 can be 
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deleted. 

8-2. Proof of Theorem 3 

It is assumed that there is a J-unitary matrix 0(k) which 
performs block triangulation as follows. 



At this time, when both sides J = (Ji©-S) - norm of the above 
expression are compared, X, Y and Z of the left side can be 
determined as follows . Where, S denote a diagonal matrix in which 
diagonal elements take 1 or -1. 
(1,1) -Block matrix 



= ^ + C t+1 (i? t+1|t -^ |t . lSP rT^ +i 

= R*t + (i^+i - Rt + i) - (12^ - i^) = r €M1 

Thus, x-^i is obtained from - J^iflSw. *« - J* 

Here, attention should be paid to the fact that 

is established. 



(2,1) -Block matrix 



X 0 
Y Z 



i _i 
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YJiX 

■ [iWVHA] 



l o 



By this , Y = ^ K ™ J H^f+i^i is obtained. 



Where, C fc = (C ifc+ i*') r 



5 



(2, 2) -Block Matrix 



= [ij^fi] +^ 1 ( i7 w^- 1 ^ T ) 

By this, 
and Z = ijt+i-RTAi is obtained. 



8-3. Proof of Theorem 4 

When an observation matrix H k has a shift characteristic and 
J = (Ji © - S), 

the following relational expression is obtained by a similar method 
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to Theorem 2. 



0 



0 



0 



Ck+\L k 



(46) 



Where, 



!?(*) = 



(47) 



Where . T,(k) T (R e>k © -R,, k ) £(*) = (R <ifc+1 © -Arj+i) 



and R r , k+1 is determined so that 2(k) T (Re.* © - R r ,k)2(k) = (Re.k+i 
© - Rr ,k+i) is established. Next, when an update expression of R r ,k+i 
is newly added to the third line of expression (46) , the following 
expression is finally obtained. 



Kk+1 

0 
0 



Where, 



0 



0 

K k 

Re* ~ C k +xL k R rJ .L k C k + l 



Ck+lLk 
Rr.k 



J 



0 



0 

Kk 



-p l >L k R- l k L T k C T k+l p- l *L k - 



0 



0 

Kk 

T »T 



(48) 



From the correspondence of the respective terms of 3 x 2 block 
matrixes of both sides, the following update expression of a gain 
10 matrix K" k is obtained. 
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0 



(49) 




(51) 
(52) 



(50) 



Industrial Applicability 

In general, in a normal civil communication equipment or 
the like, calculation is often performed at single precision in 
view of the cost and speed. Thus, as the practical state 
estimation algorithm, the present invention would have effects 
in various industrial fields. Besides, the invention can be 
applied to an echo canceller in a communication system or a sound 
system, sound field reproduction, noise control and the like. 
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CLAIMS 

1 . (Canceled) 

A system estimation method for making state estimation robust and 

optimizing — a — forgetting — coefficient — p — s imul t aneous ly — in — an 
estimation algorithm, in which 

for a state space model expres s ed by following expressions; 




10 



X fe t a state vector or simply a state, 
W fc : a system noise, 
vv: — an observation noise, 
15 Y fe t an observation signal, 

fie? dynamics of a system, an d 
G fe i a drive matrix, 

a maximum energy gain to a filter error from a disturbance 

20 weighted with — tfee — forgetting coefficient — p — as — an evaluation 
criterion i3 suppressed to be smaller than a term corresponding 
to a previously given upper limit value and 
the system estimation method comprises : 

a step at which a processing section input s the upper limit 

25 value y#, the ob s ervation signal y k a3 an input of a filter and 
a value including an observation matrix H te from a storage section 
or an input section; 

a step at which the processing s ection determines the 

forgetting coefficient p relevant to the state space model in 
30 accordance with the upper limit value y# r 

a step at which the processing section reads out an initial 

value or a value including the observation matrix Hfe-at a time from 
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the otorago section and uso9 — tho forgetting coefficient p to 
execute a hyper H co filter expressed by a following expression s 

*^ fr - m F k ^e^| k ^— <— K a ^^y k H fc F k ^x^^ k _ i -)- 

here, 

x^ k ; an estimated value of a 3tate x k at a time k U3ing observation 
-s-argnals yo to yk~r 
Kg^ k; a filter gain, 

a step at which the processing section stores an obtained 

value relating t o the hy p er filter i nto th e st orag e section; 

a — step — at — which — the — processing — section — calculates — an 

existence condition based on the upper limit value and the 
forgetting coefficient p by the obtained observation matrix H±-ar 
tho observation matrix H± and the filter gain K^, and 

a step at which the processing scotion sots the upper limit 

value to be small within a range where the existence condition 
-i-s sa tisfied at each time and stores the value into th e s torage 
section, — by decreasing the upper limit value and repeating the 
step of executing the hyper H w filter. 

2. (Amended) The system estimation method according to 
claim 4 r7 or 9 or 11 f wherein the processing section calculates 
the existence condition in accordance with a following expression: 

= S^l l + l -^-HjH i >0 1 i = 0,...,fc (17) 

3. (Amended) The system estimation method according to 
claim 7 or 9 or lli , wherein the processing section calculates 
the existence condition in accordance with a following expression: 
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here, 



-Q^i + P7/>0, i = (),•••, 



(18) 



(19) 



4. ( Amended ) The system estimation method according to 
claim 7 or 9 or 113 r, wherein the forgetting factor p and the upper 
limit value Yf have a following relation: 

0 < p = 1 - x(Yf) s 1* where x(Yf) denotes a monotonically damping 
function of Yf to satisfy xU) = 1 and x(°°) = °« 

5. (Canceled) 

10 The system estimation mothod aooording to claim 1 , whoroin at the 
stop of executing the hyper H» filter, 

the — processing — section — obtains — the — filter — gain — K 3 _ k by 

f o l lowing expression s: 



h\k = HkXk\k 

x k \ k = F k -ix k -i\k-i + K,t(yk - H k F jt_i£*-i|fc-i) 
K*,k = S k \k-iHl{H k Ek\k-iH k + p) -1 



%k\k = %k\k-l ~ ^k\k-\C k R etk CkSk\k- 

Ek+\\k = ( F kEk\kFl)/p 



here , 



e/,i = Zi\i-HiXi } xo|o = aco. £i|o--£o 

P o 



Re,k = Rk + CkEk\k-lCl> Rk - 

0 < p = 1 - xiv) < 1. 7/ > 1 
GkG T k = ^-FkSkikFl 



0 -P7/ 



C k = 



(16) 



} 



(10) 

(11) 
(12) 

(13) 



(14) 
(15) 



15 



wherein a right aide of the expression (16) — can bo mere 
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generalized, 
hero r 

the state vootor or simply tho stato, 
¥ k t the observation signal, 
g fe ; tho output signal, 
fi e- : the dynamics of the system^ 
H fe : tho observation matrix, 

x-^ fc : the estimated value of the state x k at the time k using the 
observation signals y 0 to y k -r 

Ir*^ * : corresponding to a oovarianoe matrix of an error of x" te | k -r 

K g - k ; tho filter gain , 

e f _ 3 ; the filter error, and 

Rer * s an auxiliary variable, 

6 . (Canceled) 

Th e system estimation method according to claim 5, wherein th e 
step of executing the hyper ft* filter includes : 

a step at which the processing section calculates the filter 

gain by using the expression (12) based on an initial condition; 

a step at which the processing section updates a filter 

equation of the Hc o filter of the expression (11); 

a s tep at which the proces s ing section calculates ^k j k and 

j^ifc+^l k by using the expression (13); and 

a step at which the processing s ection repeatedly execute s 

the respective steps while advancing the time k. 

7 . (Amended) The system estimation method according to 
claim 1, wherein at tho step of executing the hyper Ho> filter, 

the processing section calculates tho A system estimation 

method for making state estimation robust and optimizing a 

forgetting factor p simultaneously in an estimation algorithm, 
in which 
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for a state space model expressed by following expressions; 

x k +i = F k x k + G k w k 
Yk = H k x k + v k 
z k = H k x k 
here, 

x k ; a state vector or simply a state, 
w k ; a system noise, 
v k : an observation noise, 
y k ; an observation signal, 
10 z k : an output signal, 

F k : dynamics of a system, and 
G k ; a drive matrix, 

as an evaluation criterion, a maximum value of an energy 

gain which Indicates a ratio of a filter error to a disturbance 
15 including the system noise w k and the observation noise v k and is 
weighted with the forgetting factor p is suppressed to be smaller 
than a term corresponding to a previously given upper limit value 
If, and 

the system estimation method comprises; 

20 a step at which a processing section inputs the upper limit 

value y f , the observation signal y k as an input of a filter and 
a value including an observation matrix H k from a storage section 
or an input section; 

a step at which the processing section determines the 

25 forgetting factor p relevant to the state space model in accordance 
with the upper limit value y f ; 

a step of executing a hyper Ha, filter at which the processing 

section reads out an initial value or a value including the 
observation matrix H k at a time from the storage section and obtains 
30 a^f liter gain K Sfk by using the forgetting factor p and a gain matrix 
K k and by following expressions (20) to (22), or, expression (20) 
and expressions which are deleted Ji' 1 and Ji in the expressions 
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(21) and (22) , ; 



10 



15 



*k|* = «t-i|t-i + K t>k {y k - iffcifc-in-i) 



4 


. 1 
^t-Efcpfe-l 


0(A) = 


4* 


0 


0 


p-^k\k-l . 







(20) 
(21) 

(22) 



Where, 

Rk = RlJiRf, Rl = 



pi 0 
l 



0 pi^i 

R e j c = R k + C k E klk . l Cl, C k = 



Ji = 



1 0 
0 -1 



, 1 



, Ek\k-i = Ek\k-i S klk-i 



H k 
H k 



(23) 



0(k) denotes a J -unitary matrix, that is, satisfies 
0(k)J0H(k) T = J, J = (Ji 0 I), I denotes a unit matrix, 
K k (:,l) denotes a column vector of a first column of the 
matrix K k . 

wherein J^ T * and can be deleted in the expressions (21) 
and (22), 
here , 

x~ k | k : the estimated value of the state x k at the time k using the 

observation signals y 0 to y k , 

y k : the observation signal, 

F k : the dynamics of the system, 

K s#k : the filter gain, 

H k : the observation matrix, 

2%!*: corresponding to a covariance matrix of an error of x~ k | k , 
0(k): the J-unitary matrix, and 
R e , k : an auxiliary variable, 

a step at which the processing section stores an estimated 

value of the state x k by the hyper Hoo filter into the storage 
section; 

a step at which the processing section calculates an 
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existence condition based on the upper limit value Yf and the 
forgetting factor p by the obtained observation matrix Hj or the 
observation matrix Hj and the filter gain K s , 1# and 

a step at which the processing section sets the upper limit 

value to be small within a range where the existence condition 
is satisfied at each time and stores the value into the storage 
section, by decreasing the upper limit value Yf and repeating the 
step of executing the hyper Hop filter. 

8. (Amended) The system estimation method according to 
claim 7 , wherein the step of executing the hyper Hoo filter includes : 

a step at which the processing section calculates K k and 
2~k+i|k 1/2 by using the expression (22); 

a step at which the processing section calculates the filter 
gain K s , k based on tteean initial condition and-of 2~ k | k .i and C k , and 
the matrix gain k k by using the expressions (21) and (22) ; 

a step at which the processing section updates a filter 
equation of the Hoo filter of the expression (20); and 

a step at which the processing section repeatedly executes 
the respective s tops the step of calculating by using the expression 
(20), the step of calculating by using the expressions (21) and 
(22), and, the step of updating while advancing the time k. 

9 . (Amended) The 3y3tem estimation m et hod according to 
claim 1, wherein at the stop of executing the hyper Hoo filter, 

the processing s ection A system estimation method for making 

state estimation robust and optimizing a forgetting factor p 
simultaneously in an estimation algorithm, in which 

for a state space model expressed by following expressions: 

x k + x = F k x k + G k w k 
y k = H k x k + v k 
z k ■ H k x k 
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here, 

x k : a state vector or simply a state, 

w k : a system noise, 

v k : an observation noise, 

y k ; an observation signal, 

z k ; an output signal, 

F k ; dynamics of a system, and 

G k : a drive matrix, 

as an evaluation criterion, a maximum value of an energy 

gain which indicates a ratio of a filter error to a disturbance 
including the system noise w k and the observation noise v k and is 
weighted with the forgetting factor p is suppressed to be smaller 
than a term corresponding to a previously given upper limit value 
Yf * and 

the system estimation method comprises: 

a step at which a processing section inputs the upper limit 

value Yf/ the observation signal y k as an input of a filter and 
a value including an observation matrix H k from a storage section 
or an input section; 

a step at which the processing section determines the 

forgetting factor p relevant to the state space model in accordance 
with the upper limit value y±; 

a step of executing a hyper H«> filter at which the processing 

section reads out an initial value or a value including the 
observation matrix H k at a time from the storage section and obtains 
the a filter gain K s , k by using the forgetting factor p and a gain 
matrix K k and by following expressions: 
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«*|* = *t-i|t-i + K,jk{y k - JTfcXfc^it-i) 
K $M = l)/i^(l, 1) , K k = ^(J^Rjj)^ 

I 



(61) 
(62) 

0(Jk) (63) 



here, 0(k) denotes an arbitrary J-unitary matrix, and C k = 
Cv + iVP is established, where 

R k = R\j l R][, R\ = 



pi o 1 fiol - £ 

0 J ' Jl = [ 0 -1 J ' = ^W^V* 

fl e ^ = H Jb + C 4 i7 4|Jb _ 1 C^ l C * = [ ^* ] > ^ = R^R^, x 0 |o = x 0 (23) 



here, 

x~ k |k: the estimated value of the state x k at the time k using the 
observation signals y 0 to y k , 
y k : the observation signal, 
K S/k : the filter gain, 
H k : the observation matrix, 
0(k): the J-unitary matrix, and 
Re, k : an auxiliary variable. 
10 a step at which the processing section stores an estimated 

value of the state x k by the hyper Hoo filter into the storage 
section; 

a step at which the processing section calculates an 

existence condition based on the upper limit value Yf and the 
15 forgetting factor p by the obtained observation matrix Hj or the 
observation matrix Hi and the filter gain K s<1 , and 

a step at which the processing section sets the upper limit 

value to be small within a range where the existence condition 
is satisfied at each time and stores the value into the storage 
20 section, by decreasing the upper limit value Yf and repeating the 
step of executing the hyper Hoo filter. 
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10. (Amended) The system estimation method according to 
claim 9 , wherein the step of executing the hyper Hoo filter includes : 

a step at which the processing section calculates K~ k based 
on an initial condition of R e ,k+i, Rr.k+i and L\+ lf and the matrix 
gain K\ by using the expression (63); 

a step at which the processing section calculates the filter 
gain K^^-K s , k based on the initial condition and by using the 
expression (62); 

a step at which the processing section updates a filter 
equation of the Hoo filter of the expression (61); and 

a step at which the processing section repeatedly executes 
— respective — s^ep ethe step of calculating by using the 
expressions (62) and (63), the step of calculating by using the 
expression (61) , and, the step of updating while advancing the 
time k. 

11 • (Amended) The oyotom estimation method according to 
claim 1, wherein at the step of executing the hyper R» filter, 

the processing section A system estimation method for making 

state estimation robust and optimizing a forgetting factor p 
simultaneously in an estimation algorithm, in which 

for a state space model expressed by following expressions: 

x k +i = F k x k + Gfc Wk 
Yk = H k x k + v k 
Zk = H k x k 
here, 

x k ; a state vector or simply a state, 

w k : a system noise, 

v k : an observation noise, 

y k : an observation signal, 

z k : an output signal, 

F k : dynamics of a system, and 
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Gk: a drive matrix, 

as an evaluation criterion, a maximum value of an energy 

gain which indicates a ratio of a filter error to a disturbance 
including the system noise w k and the observation noise v k and is 
weighted with the forgetting factor p is suppressed to be smaller 
than a term corresponding to a previously given upper limit value 
If . and 

the system estimation method comprises: 

a step at which a processing section Inputs the upper limit 

value Yf, the observation signal y k as an input of a filter and 
a value including an observation matrix H k from a storage section 
or an input section; 

a step at which the processing section determines the 

forgetting factor p relevant to the state space model in accordance 
with the upper limit value Y f; 

a step of executing a hyper Ho> filter at which the processing 

section reads out an initial value or a value Including the 
observation matrix H k at a time from the storage section and obtains 
the a filter gain K s , k by using the forgetting factor p and a gain 
matrix K" k and by following expressions: 
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- P -n k R-iL T k c T k +1 



Kk+i 




0 


0 







Ljfc+l = p il/jt - 



0 



-R e ,*+1 = R e,k ~ Ck+lLkRrjkLkC^ 



-t v r 



Where, 



(25) 
(26) 

(27) 

(28) 

(29) 
(30) 



i+1 



fte.l = H] + C^noCf, ft, = 

e ft,^ 



1 0 
0 0 
0 1 



, H k+1 = [u k+l u(k + l-N)} = [u(k + l)u k ], #, = [«(1),0,...,0] 

, A|o = diag{ /9 2 ,p 3 ,...,/+2j ) p= i- x ( 7/ ) 

#o = 0, i 0 |o = x 0> K t = p"*K t (31) 



P 0 



0 -r* 



-1 0 

0 p~" 



wherein tho above expressions oan bo arranged alao with 
roopodt to Kje inotoad of K" H - 
here, 

y k : the observation signal, 
F k : the dynamics of the system, 
H k : the observation matrix, 

x~ k |k: the estimated value of the state x k at the time k using the 
observation signals y 0 to y k , 

K s , k : the filter gain, obtained from the gain matrix K" k , and 
Re.k, L~k: an auxiliary variable. 

a step at which the processing section stores an estimated 

value of the state x k by the hyper Ho, filter into the storage 
section; 

a step at which the processing section calculates an 

existence condition based on the upper limit value y f and the 
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forgetting factor p by the obtained observation matrix Hj or the 
observation matrix Hi and the filter gain K s , ± , and 

a step at which the processing section sets the upper limit 

value to be small within a range where the existence condition 
is satisfied at each time and stores the value into the storage 
section, by decreasing the upper limit value y f and repeating the 
step of executing the hyper H«, filter. 

12. (Canceled) 

The system estimation method according to claim 11, wherein the 
step of executing the hyper II «> filter includes; 

a — step — at — which — the — processing — section — recursively 

calculates K~W based on a previously determined initial condition 
and by using the expression (27); 

a step at which the processing section calculates the system 

gain K STk by using the expression (26 ); 

a step — at — which — the proces -sing — section — calculates — the 

existence oonditlon; - 

a step at which the processing section update s a filter 

equation of the IIo» filter of the expression (25) when the existence 
condition is satisfied, and repeatedly executes the respective 
stops repeatedly while advancing the time k; and 

a step of increasing the upper limit value y± when the 

existence condition is not satisfied. 

13. (Amended) The system estimation method according to 
claim - 17 or 9 or 11 , wherein an estimated value z v k | k of the output 
signal is obtained from the state estimated value x~ k | k at the time 
k by a following expression: 

z v k | k = H k x~ k | k . 

14. (Amended) The system estimation method according to 
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claim 7 or 9 or wherein the Hoo filter equation is applied 

to obtain the state estimated value x~ k | k =[h~i[k] , • " , h~ N [k] ] T 
a pseudo-echo is estimated by a following expression: 

i = E * = 0,1,2,... (34) 

and 

an echo canceller is realized by canceling an actual echo 
by the obtained pseudo-echo. 

15. (Amended) A system estimation program for causing a 
computer to make state estimation robust and to optimize a 

forgetting factor p simultaneously in an estimation algorithm, 
in which 

for a state space model expressed by following expressions : 

x k +i = F k x k + G k Wk 
y k = H k x k + v k 
z k = H k x k 
here, 

x k : a state vector or simply a state, 

w k ; a system noise, 

v k ; an observation noise, 

y k : an observation signal, 

z k ; an output signal, 

F k ; dynamics of a system, and 

G k ; a drive matrix, 

as an evaluation criterion, a maximum value of an energy 

gain which Indicates a ratio of a filter error to a disturbance 
including the system noise w k and the observation noise v k and is 
weighted with the forgetting factor p is suppressed to be smaller 
than a term corresponding to a previously given upper limit value 
Xt , and 
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for a state space model expressed by following expression s t 



here, 

x k ; a state vector or simply a state, 

Wfci a system noise, 

v ^; an observation noise, 

Yk: an observation signal, 

-&k: an output signal, 

fie . : dynamics of a system, and 

Gict a drive matrix^ 

a maximum energy gain to a filter error from a disturbance 

weighted with the forgetting coefficient p as an evaluation 
criterion is suppressed to bo smaller than a term corresponding 
to a previously given upper limit value y^ r and 

the system estimation program causes the computer to 
execute : 

a step at which a processing section inputs the upper limit 

value Vf, the observation signal y k as an input of a filter and 
a value including an observation matrix H k from a storage section 
or an input section; 

a step at which the processing section determines the 

forgetting factor p relevant to the state space model in accordance 
with the upper limit value Y f? 

a step of executing a hyper Hop filter at which the processing 

section reads out an initial value or a value including the 
observation matrix H k at a time from the storage section and obtains 

a filter gain K s , k by using the forgetting factor p and a gain matrix 
K" k and by following expressions; 
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%k\k = &k-\\k-i + K 3ik (y k - if jbX*-i|jt_i) 
K,, k = piK k (:A)/IUAhl) 







0 


0 







-p *L k R r lL k C M 



Lk+i = P~*L k - 



0 



>-i- t ~ t 



Where , 



(25) 
(26) 

(27) 

(28) 

(29) 
(30) 



J 



L 0 = 



1 0 
0 0 
0 1 



. ^Hi = k+itt(A + l-iV)] = [u(i + l)tt t ], Hi=[u(l) ) 0,...,0] 

, i7 1|0 = diag{p 2 ,p 3 ) ...,p^ }> p= i_ x ( T/ ) 

, #0 = 0, x 0 |o = x 0 , Kt=p- { K k (31) 



6 ji , jb = 



P 0 
0 -pyj 



-1 0 

0 p~ N 



here, 

y k : the observation signal, 

F k ; the dynamics of the system, 

H k ; the observation matrix, 

x" k | k ; the estimated value of the state x k at the time k using the 
observation signals y 0 to y k , 

K s , k ; the filter gain, obtained from the gain matrix K" k , and 
Re.k> L\: an auxiliary variable. 

a step at which the processing section stores an estimated 

value of the state x k by the hyper H, filter into the storage 
section; 

a step at which the processing section calculates an 

existence condition based on the upper limit value V f and the 
forgetting factor p by the obtained observation matrix Hi or the 
observation matrix H ± and the filter gain K s ,i, and 
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a step at which the processing section sets the upper limit 

value to be small within a range where the existence condition 
is satisfied at each time and stores the value into the storage 
section, by decreasing the upper limit value Yf and repeating the 
step of executing the hyper Hoo filter. 

a stop at which a proc e ssing soction inputs the uppor limit 

value Yz, the observation signal y k as an input of a filter and 
a value including an observation matrix H k from a storage section 
or an input section; 

a step — at — which — the processing — section — determines — the 

forgetting coefficient p relevant to the state s pace model in 
accordance with the uppor limit value y& r 

a stop at which the processing section roads out an initial 

value or a value including the observation matrix H^at a time from 
the storage section and U3cs the forgetting coefficient p to 
execute a hyper Ho o filter expressed by a following expression ; 

i HfcFfc-^x^x l[k in - 
here, 

3^ic| k .; an estimated value of a state x k at a time k using observation 

signals y 0 to y k -r 

fi t ; dynamics of the system, and 

K gTk ; a filter gain> 

a step at which the processing section stores an obtained 

value relating to the hyper KUo filter into th e storage section; 

a — step — at — which — the — processing — section — calculates — an 

existence condition ba s ed on the upper limit value y* and th o 
forgetting coefficient p by the obtained observation matrix H± — or 
the observation matrix H± and the filter gain K^-r-^and 

a stop at which the processing section sets the upper limit 

value to be small within a range where the existence condition 
- is satisfied at each time and stores the value into the storage 
section by decreasing the uppor limit value and repeating the 
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stop of Qxocuting the hyper Ho > filter. 

16. (Amended) A computer readable recording medium 
recording a system estimation program for causing a computer to 
make state estimation robust and to optimize a forgetting factor 
p simultaneously in an estimation algorithm, in which 

for a state space model expressed by following expressions; 

Xk+i = *kX k + G k w k 
y k = H k x k + v k 
z k = H k x k 
here, 

x k ; a state vector or simply a state, 

w k ; a system noise, 

v k : an observation noise, 

y k : an observation signal, 

z k : an output signal, 

F k ; dynamics of a system, and 

G k : a drive matrix, 

as an evaluation criterion, a maximum value of an energy 

gain which indicates a ratio of a filter error to a disturbance 
including the system noise w k and the observation noise v k and is 
weighted with the forgetting factor p is suppressed to be smaller 
than a term corresponding to a previously given upper limit value 
I*, and 

for a state 3pace model expressed by following expression s : 

Xfc+± - i G k w te 

Yk a HfcXjt i v k 

here, 

X fe t a state vector or simply a state, 

Wfe-t — a system noise, 

V k ; an observation noise, 
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y k : an observation signal, 
^ k : an output signal, 
F fe t dynamics of a system, and 
G k ; a drive matrix, 

a maximum energy gain to a filter error from a disturbance 

weighted with the forgetting coefficient p as an evaluation 
criterion i s suppressed to be smaller than a term corre s ponding 
to a previously given upper limit value v* t— etn4 

the computer readable recording medium recording the system 
estimation program causes the computer to execute: 

a step at which a processing section inputs the upper limit 

value Vf, the observation signal y k as an input of a filter and 
a value including an observation matrix H k from a storage section 
or an input section; 

a step at which the processing section determines the 

forgetting factor p relevant to the state space model in accordance 
with the upper limit value Vf: 

a step of executing a hyper Hop filter at which the processing 

section reads out an initial value or a value including the 
observation matrix H k at a time from the storage section and obtains 
a filter gain K s , k by using the forgetting factor p and a gain matrix 
K" k and by following expressions; 
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K, <k = p*K k (:,l)/R e +(l,l) 

-p-sl k R-lL T k C T k+l 



K k +\ 




0 


0 




.K k . 



r -' - 



0 



- -Rr,fc ~ L k C k+l R ek C k+i L k 

Where, 



(25) 
(26) 

(27) 

(28) 

(29) 
(30) 



Hk+i 



L 0 = 



1 0 
0 0 
0 1 



. &k + i = [u M u(k + l-N)] = [u(k + l)u k ) t H 1 = [u(l),0,...,0] 

, ^i|o = diag{p 2 , /9 3 ,...,p w +2} ) p=l- x (y f ) 

Ko = 0, x 0 |o = x 0 , = (31) 



Re,i = Ri 4 C,27i| 0 C, , ft, = 

6 ^ o = 



P 0 
0 



-1 0 

0 p~ N 



here, 



y k ; the observation signal, 

F k ; the dynamics of the system, 

H k ; the observation matrix, 

x" k | k ; the estimated value of the state x k at the time k using the 
observation signals y 0 to y k , 

K s , k ; the filter gain, obtained from the gain matrix K' k , and 
Re.k, L"k: an auxiliary variable. 

a step at which the processing section stores an estimated 

value of the state x k by the hyper H„ filter into the storage 
section; 



a step at which the processing section calculates an 



existence condition based on the upper limit value Vf and the 
forgetting factor p by the obtained observation matrix Hj or the 
observation matrix Hj and the filter gain K s , lt and 
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a step at which the processing section sets the upper limit 

value to be small within a range where the existence condition 
is satisfied at each time and stores the value into the storage 
section, by decreasing the upper limit value Vf and repeating the 

step of executing the hyper Hg> filter. a — &tep — at — which — a 

processing section inputs the upper limit value y#, the observation 
signal y^— -as — an input of a filter and a value including an 
observation matrix H te - from a storage section or an input section; 

a 3 top — at which the proce s sing section determines — the 

forgetting coefficient p relevant to the state space model in 
accordance with the upper limit value y^ -r 

a step at which the processing section roads out an initial 

value or a value including the observation matrix H fe at a time from 
the storage section and uses the forgetting coefficient p to 
execute a hyper H « > filter expressed by a following expression ; 
x^ fo m Fte ^x^^lfe^— ^^MYh HfeFfc-^ x afr in - 
here # 

an estimated value of a state at a time k using observation 
signals y& to y^ — 
fi e ; dynamics of the system, and 
K* ^; a filter gain, 

a step at which the processing section stores an obtained 

value relating to the hyper Hg > filter into the storage section; 

a — step — a* — which — the — processing — section — calculates — an 

existence condition based on the upper limit value y € and the 
forgetting coefficient p by the obtained observation matrix Ki-er 
the observation matrix H± and the filter gain K^,^— and 

a stop at which the processing section sets the upper limit 

value to bo 3mall within a range where the existence condition 
is satisfied at each time and stores the value into the storage 
section by decreasing the upper limit value y* and repeating the 
step of executing the hyper H<* filter . 
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17. (Amended) A system estimation device for making state 
estimation robust and optimizing a forgetting factor p 
simultaneously in an estimation algorithm, in which 

for a state space model expressed by following expressions: 

x k +i = F k x k + G k w k 
y k = H k x k + v k 
Zk = H k x k 
here, 

x k ; a state vector or simply a state, 

w k ; a system noise, 

v k : an observation noise, 

y k : an observation signal, 

z k : an output signal, 

F k : dynamics of a system, and 

G k ; a drive matrix, 

as an evaluation criterion, a maximum value of an energy 

gain which indicates a ratio of a filter error to a disturbance 
including the system noise w k and the observation noise v k and is 
weighted with the forgetting factor p is suppressed to be smaller 
than a term corresponding to a previously given upper limit value 
It . and 

for a state s pace model expressed by following expressions; 

Xfe+i ■ FfcXfc — 1 — GfcWfc 

hero, 

Xfc-: — a state vector or simply a state, 

a system noise, 
Vh: an observation noise, 
y ^ ; an observation signal, 
^ k t an output signal. 
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fi c t dynamics of a system -; — a**d 
6k: a drive matrix, 

a maximum energy gain to a filter error from a disturbance 

weighted with the forgetting coefficient p as an evaluation 
criterion io suppressed to be smaller than a term corresponding 
to a previously given upper limit value y #t 

the system estimation device comprises: 

a processing section to execute the estimation algorithm; 

and 

a storage section to which reading and/or writing is 
performed by the processing section and which stores respective 
observed values, set values, and estimated values relevant to the 
state space model, 
further comprising: 

a means at which the processing section inputs the upper 

limit value Vf , the observation signal y k as an input of a filter 
and a value including an observation matrix H k from the storage 
section or an input section; 

a means at which the processing section determines the 

forgetting factor p relevant to the state space model in accordance 
with the upper limit value y f ; 

a means of executing a hyper Hoo filter at which the processing 

section reads out an initial value or a value including the 
observation matrix H k at a time from the storage section and obtains 
a filter gain K s , k by using the forgetting factor p and a gain matrix 
K" k and by following expressions: 
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*k\k ~ *k-\\k-\ + KaJk{Vk ~ HkZk-l\k-l) 

K Stk = p^(:,l)/^(l,l) 

-p^L k R;yL T k C T k+1 



Kk+l 




0 


0 




.*k. 



L k +\ = p *L k 



0 

[Kk\ 



1 w " 

R e,k C k+lL k 



u - _j - T w T 
-Re.k+1 = -Re,fc - Ck+\LkRr >k L k C k+l 

#r,fc+l = Rrjt ~ L k G k+l R e,k C k+l L k 



(25) 
(26) 

(27) 

(28) 

(29) 
(30) 



Where , 



J 



1 0 
0 0 
0 1 



. = K+i «(* + 1 - AO] = [u(k + 1) u*], Hi = («(1), 0 0] 

27,(0 = diag{p 2 , />',..., p = 1 - X ( 7/ ) 

, #0 = 0, x 0|0 = i„, Ki = p' k K k (31) 



6 ^ 



P o 

PI? J 



0 -~* 2 



-1 0 

o P -»\ 



here, 

y k ; the observation signal, 

F k : the dynamics of the system, 

H k ; the observation matrix, 

x~ k | k ; the estimated value of the state x k at the time k using the 
observation signals y 0 to y k , 

K Sjk : the filter gain, obtained from the gain matrix K" k , and 
Re.R* L" k : an auxiliary variable, 

a means at which the processing section stores an estimated 

value of the state x k by the hyper Ho filter into the storage 
section; 

a means at which the processing section calculates an 

existence condition based on the upper limit value Vf and the 
forgetting factor p by the obtained observation matrix Hi or the 
observation matrix Hi and the filter gain K s .j, and 
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a means at which the processing section sets the upper limit 

value to be small within a range where the existence condition 
is satisfied at each time and stores the value into the storage 
section, by decreasing the upper limit value Vf and repeating the 
means of executing the hyper Hoo filter, — the proces s ing section 
inputs the upper limit value y ^ , the observation signal y fe a s an 
input of a filter and a value including an observation matrix ¥L& 
from a storage section or an input section , 

the proce s sing section determi nes the forgetting 

10 coefficient p relevant to the state space model in accordance with 
the upper limit value y# t 

the processing section reads out an initial value or a value 

including the observation matrix H k at a time from the storage 
section and use s the forgetting coefficient p to execute a hyper 
15 Hoo filter expressed by a following expression : 

here, 

x^-fefc; an estimated value of a state x k at a time k using ob s ervation 
signals y< > to y te - 
20 fi e : dynamics of the system, — and 
a filter gain, 

the proce s sing section store s an obtained value relating 

to the hyper Hoo filter into the storage section, 

the processing s ection calculates an existence condition 

25 based on the upper limit value and the forgetting coefficient 
p by the obtained observation matrix H± or the observation matrix 
H± and the filter gain K g-± , and 

the processing part sets the upper limit value to be small 

within a range where the existence condition is satisfied at each 
30 time and stores the value into the storage s ection by decreasing 
the upper limit value repeating the 3tcp of executing the 

hyper H oo filter. 
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